/*
Copyright (c) 2008 Chagrin Research, LLC

Permission is hereby granted, free of charge, to any person
obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without
restriction, including without limitation the rights to use,
copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following
conditions:

The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
OTHER DEALINGS IN THE SOFTWARE.
*/

#include <assert.h>
#include <stddef.h>
#include <math.h>
#include <arpa/inet.h>
#include "primitives.h"
#include <unistd.h>
#include "utm.h"

#define NRS(x) (x)=ntohs((x))
#define NRL(x) (x)=ntohl((x))

const double meters_per_nm = 1852.0;

typedef enum special_code
{
NUMBER=-1,
BLANK=0,
TH,
ND,
RF
} special_code;

double get_engineering_units(unsigned short code, special_code* sc)
{
	double res = 0.0;
	*sc = NUMBER;
	if(0x8000 & code)
	{
		switch(((unsigned char)code))
		{
			case 0x00:
				*sc = BLANK;
				break;
			case 0x01:
				*sc = TH;
				break;
			case 0x02:
				*sc = ND;
				break;
			case 0x03:
				*sc = RF;
				break;
			default:
				assert(0);
		};
	}
	else
	{
		res = (double)((char)code);
		if(0x4000 & code)
		{
			//data in least significant byte is scaled by 100
			res /= 100.0;
		}
		else
		{
			if(0x2000 & code)
			{
				//data in least significant byte is scaled by 20
				res /= 20.0;
			}
			else
			{
				if(0x1000 & code)
				{
					//data in least significant byte is scaled by 10
					res /= 10.0;
				}	
			}
		}
		
		if(0x0100 & code)
			res = -1.0 * res;
		
	}
	return res;
}

//range from 5 dBZ to 75 dBZ
rgba8 levels[] = {
{66,66,66,255},{99,99,99,255},
{40,126,40,255},{60,160,20,255},{120,220,20,255},
{250,250,20,255},{250,204,20,255},{250,153,20,255},
{250,79,20,255},{250,0,20,255},{220,30,70,255},{200,30,100,255},
{170,130,50,255},{255,0,156,255},{255,255,255,255}};

void translate_levels(int width,int i,double* grid_data,rgba8* row)
{
	int j;
	for(j=0;j<width;j++)
	{
		double value = grid_data[j];
		if(value > 75.0)
			value=75.0;
		else
		{
			if(value < 0.0)
				value=0.0;
		}
				
		value = 16.0 * (value/75.0);
		row[j] = levels[(int)value];
		row[j].a = 128;
		
		if(value==0.0)
			row[j].a = 0;
	}
}

//Level III data with raster data format uses FlatEarth Projection. This projection surface is tangent at some point (lat0, lon0) and has a y axis rotated from true North by some angle. We call it "flat" because it should only be used where the spherical geometry of the earth is not significan
static grid* process_raster_data(void* data, msg_header* m,double meters_per_pix)
{
	NRS(m->raster_data.packet_code2);
	NRS(m->raster_data.packet_code3);
	NRS(m->raster_data.i_coordinate_start);
	NRS(m->raster_data.j_coordinate_start);
	NRS(m->raster_data.x_scale_int);
	NRS(m->raster_data.x_scale_fractional);
	NRS(m->raster_data.y_scale_int);
	NRS(m->raster_data.y_scale_fractional);
	NRS(m->raster_data.number_of_rows);
	NRS(m->raster_data.packing_descriptor);
	
	grid* g = grid_create(m->raster_data.number_of_rows,m->raster_data.number_of_rows);
	
	double station_lat = m->latitude*0.001;
	double station_lng = m->longitude*0.001;
	
	double station_x,station_y;
	
	LLtoUTM(WGS_84,station_lat,station_lng,&station_y,&station_x,g->zone);
	
	double dlat,dlng;
	double lat,lng;
	
	UTMtoLL(WGS_84,station_y - 1000.0,station_x,g->zone,&lat,&lng);
	
	dlat = station_lat - lat;
	
	UTMtoLL(WGS_84,station_y,station_x - 1000.0,g->zone,&lat,&lng);
	
	dlng = station_lng - lng;
	
	assert(dlat > 0.0 && dlng > 0.0);
	
	double degrees_per_meter_lat = dlat / 1000.0;
	double degrees_per_meter_lng = dlng / 1000.0;
	
	double width = (double)m->raster_data.number_of_rows * meters_per_pix;//width in meters
	
	lat = station_lat - (width/2.0) * degrees_per_meter_lat;
	lng = station_lng - (width/2.0) * degrees_per_meter_lng;
	
	grid_set_geometry(g,lat,lng,meters_per_pix,meters_per_pix);
	
	g->bounds.southwest_lat = lat;
	g->bounds.southwest_lng = lng;
	
	lat = station_lat + (width/2.0) * degrees_per_meter_lat;
	lng = station_lng + (width/2.0) * degrees_per_meter_lng;
	
	g->bounds.northeast_lat = lat;
	g->bounds.northeast_lng = lng;
	
	data = data + sizeof(raster);
	
	int i,j;
	for(i=0;i<m->raster_data.number_of_rows;i++)
	{
		unsigned short bytes_in_line;
		memcpy(&bytes_in_line,data,sizeof(bytes_in_line));
		NRS(bytes_in_line);
	
		data += sizeof(bytes_in_line);
		int ii;
		j = 0;
		for(ii=0;ii < bytes_in_line/2;ii++)
		{
			unsigned short t;
			memcpy(&t,data,sizeof(unsigned short));
			NRS(t);
			rle tmp;
			memcpy(&tmp,&t,sizeof(rle));
			data += sizeof(unsigned short);
			
			int iii;
			
			unsigned short value0 = m->data_level_threshold[tmp.color0];
			unsigned short value1 = m->data_level_threshold[tmp.color1];
			special_code sc = BLANK;
			double eu0 = get_engineering_units(value0,&sc);
			double eu1 = get_engineering_units(value1,&sc);
			
			if(eu0 > g->max_value)
				g->max_value = eu0;
			if(eu1 > g->max_value)
				g->max_value = eu1;
			
			for(iii=0;iii<tmp.run0;iii++)
			{
				grid_set_pixel(g,j,i,eu0);
				j++;
			}
			for(iii=0;iii<tmp.run1;iii++)
			{
				grid_set_pixel(g,j,i,eu1);
				j++;
			}
		}
	}
	
	return g;
}

const char echo_tops_url[] = "curl ftp://tgftp.nws.noaa.gov/SL.us008001/DF.of/DC.radar/DS.p41et/SI.%s/sn.last -o cr.%s.last";
const char composite_url[] = "curl ftp://tgftp.nws.noaa.gov/SL.us008001/DF.of/DC.radar/DS.p37cr/SI.%s/sn.last -o cr.%s.last";

grid* nexrad(char* station, double* station_lat,double* station_lng,bounding_box* b, char* url)
{
	mapped_file* mf;
	char* header;
	double meters_per_pix=0.0;	
	char station_cmd[8192];
	
	snprintf(station_cmd,8192,composite_url,station,station);
	meters_per_pix = 0.54 * meters_per_nm;

	if(system(station_cmd))
	{
		printf("No such station");
		exit(-1);
	}
	
	snprintf(station_cmd,8192,"cr.%s.last",station);
	
	mf=mapfile(station_cmd,0);
	header = strstr(mf->data,"\r\r\n");
	header ++;
	header = strstr(header,"\r\r\n");
	
	void* data = header + 3;
	
	msg_header* m = malloc(sizeof(msg_header));
	
	memcpy(m,data,sizeof(msg_header));
	
	NRS(m->message_code);
	NRS(m->julian_date_code);
	NRL(m->timestamp);
	NRL(m->length);
	NRS(m->source_id);
	NRS(m->dest_id);
	NRS(m->number_blocks);
	
	NRS(m->divider);
	NRL(m->latitude);
	NRL(m->longitude);
	NRS(m->height);
	NRS(m->product_code);
	NRS(m->operational_mode);
	NRS(m->volume_coverage_pattern);
	NRS(m->sequence_number);
	NRS(m->volume_scan_number);
	NRS(m->volume_scan_date);
	NRL(m->volume_scan_start_time);
	NRS(m->product_gen_date);
	NRL(m->product_gen_time);
	NRS(m->elevation_angle);
	int i;
	for(i=0;i<23;i++)
		NRS(m->data_level_threshold[i]);
	
	NRL(m->offset_to_symbology);
	NRL(m->offset_to_graphic);
	NRL(m->offset_to_tabular);
	
	NRS(m->divider2);
	NRS(m->block_id);
	NRL(m->block_length);
	NRS(m->number_of_layers);
	NRS(m->layer_divider);
	NRL(m->length_data_layer);
	NRS(m->packet_code);
	
	grid* g = NULL;
	
	switch(m->packet_code)
	{
		default:
			assert(0);
			break;
		case 0xAF1F:
			//process_radial_data(data + offsetof(msg_header,radial_data),m);
			break;
		case 0xBA07:
			g = process_raster_data(data + offsetof(msg_header,raster_data),m,meters_per_pix);
			break;
			
	};
	
	*station_lat = (double)m->latitude * 0.001;
	*station_lng = (double)m->longitude * 0.001;
	
	double northing;
	double easting;
	LLtoUTM(WGS_72,*station_lat,*station_lng,&northing,&easting,g->zone);
	
	UTMtoLL(WGS_72,northing,easting,g->zone, station_lat, station_lng);
	
	b->south_west.lat = g->bounds.southwest_lat;
	b->south_west.lng = g->bounds.southwest_lng;
	b->north_east.lat = g->bounds.northeast_lat;
	b->north_east.lng = g->bounds.northeast_lng;
	
	free(m);
	unmapfile(mf);
	
	return g;
}


//ftp://tgftp.nws.noaa.gov/SL.us008001/DF.of/DC.radar/DS.p19r0/SI.kcle
grid* get_nexrad_data(char* station, double* station_lat, double* station_lng,bounding_box* b,char* url)
{
	return nexrad(station,station_lat,station_lng,b,url);
}

////////////////////

mapped_file* mapfile(const char* filename, size_t length)
{
	mapped_file* mf = malloc(sizeof(mapped_file));
	memset(mf,0,sizeof(mapped_file));
	
	mf->fd = open(filename,O_RDONLY);
	if(mf->fd > 0)
	{
		fstat(mf->fd,&mf->stat);
	
		if(!length || (mf->stat.st_size < length))
			length = mf->stat.st_size;
		
		mf->data = mmap(NULL,length,PROT_READ,MAP_PRIVATE,mf->fd,0);
		if(mf->data == MAP_FAILED)
		{
			close(mf->fd);
			free(mf);
			return NULL;
		}
		mf->length = length;
		mf->offset=0;
	}
	else
	{
		free(mf);
		return NULL;
	}
	return mf;
}

void unmapfile(mapped_file* f)
{
	munmap(f->data,f->length);
	close(f->fd);
	free(f);
}
